Families of solitons in Bragg supergratings 
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We study fundamental optical gap solitons in the model of a fiber Bragg grating (BG), which is subjected to 
i— i a periodic modulation of the local reflectivity giving rise to a supergrating. In addition, the local refractive 

00 ... . 

p ^ index is also periodically modulated with the same period. It is known that the supergrating opens an infinite 

system of new bandgaps in the BG's spectrum. We use a combination of analytical and computational 
methods to show that each emerging bandgap is filled with gap solitons (GSs), including asymmetric ones 
and bound states of the GSs. In particular, bifurcations of the GSs created by the supergrating are studied 
in terms of a geometric analysis. 
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1. Introduction 

Bragg gratings (BGs) are light-controlling structures produced by a periodic variation of the refractive 
index along an optical fiber or waveguide. Devices based on fiber gratings, such as dispersion compensators, 
sensors, elements of laser cavities, etc., are widely used in optical systems [2T] . Gap solitons (GSs) in fiber 
gratings, alias BG solitons [TTJ]. are supported through the balance between the BG-induced dispersion, 
which incorporates the bandgap in the system's linear spectrum, and the Kerr nonlinearity of the fiber or 
waveguide. Analytical solutions for BG solitons in the standard fiber-grating model are well known [HE]. 
Slightly more than half of the analytical found family is stable, as first demonstrated, within the framework 
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of the variational approximation, in Ref. |26| . and then in a consistent numerical form in Ref. Later, 
an accurate numerical technique for simulations of GSs was elaborated in Ref. [9], and a rigorous analysis 
of their stability was developed in Ref. [5] 

Following the theoretical prediction, BG solitons were created in the experiment, using a relatively short 
BG written in the cladding of an optical fiber [15] and a virtual optically-induced grating |39j . In particular, 
an essential physical achievement was the creation of slow BG solitons moving at the velocity much lower 
than the velocity of light c in vacuum (namely, 0.16c |27|). 

GSs are one of fundamental species of solitons in optics, as well as in other nonlinear media. In particular, 
the GS concept was extended to Bose- Einstein condensates (BECs), where such solitons were predicted in 
a condensate loaded into a periodic potential (optical lattice) [2] [U [3T] , and then created in the experiment 

An issue of great significance to the fundamental studies and applications is the development of methods 
for control of BG solitons. One of them is apodization [TU [25l [27J [36], i.e., the use of a grating with the 
Bragg reflectivity gradually varying along the fiber (or waveguide). In particular, a possibility was predicted 
to slow down the soliton and eventually bring it to a halt in a properly apodized BG |25j . Experimentally, it 
has been shown that the apodization helps to couple solitons into fiber gratings |14j , and may also split them 
|36j . The retardation of the BG solitons to the above-mentioned low velocity, 0.16c, was also performed 
with the help of the apodization [57] . 

Technologies which make it possible to fabricate regular fiber gratings with periodic apodization, i.e., 
an effective superlattice built on top of the BG, are well established |35j . Analysis of the light propagation 
in the so developed supergratings was performed in [5]. It was shown that the supergrating gives rise to 
extra gaps in the system's spectrum ("Rowland ghost gaps"). Solitons in the gaps were looked for in an 
approximation based on "supercoupled-mode equations" , which, essentially, assume that a soliton amounts 
to a slowly varying envelope of the supergrating's Bloch function(s). As is known, such an approximation 
(the averaging method) generally applies to the description of GSs near bandgap edges [I] . A related problem 
was considered in [2HIHI], viz., GSs in the BEC loaded into a double-periodic optical lattice. The period 
doubling opens an additional narrow "mini- gap", where stable GSs exist. 

The objective of our work is to investigate the existence of GSs in harmonic superlattices created on top 
of the usual BG in fibers or waveguides with the intrinsic Kerr nonlinearity, extending the analysis initiated 
in [49] . Our model has periodic potentials in one spatial dimension. Analytical approaches for GSs in two 
or more spatial dimensions were also developed recently in [121 G31 023 132] ■ 

The paper is organized as follows. We first give our mathematical model of the supergrating and analyze 
the bandgap structure in the model, using its dispersion relation, in Section[2] This constitutes a non-trivial 
extension of the classical spectral theory for the Mathieu equation [37], and predicts parameter values at 
which GSs may be expected to exist in the model. The sequence of new bandgaps not present in the 
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unmodulated BG is produced. 

In Sections [3] and [4] we establish the existence of GSs in the newly opened bandgaps, using a combination 
of analytical and numerical techniques. We find various families of GSs in the supcrgrating, including 
those which do not exist in the ordinary GS model (in particular, stable asymmetric solitons, and multi- 
soliton bound states). Using the Melnikov method, we first prove that GSs in the original BG model have 
counterparts in the extended model. A similar treatment was used for a perturbed Hamiltonian partial 
differential equations in [^Hj. Then, the method of averaging is used to demonstrate the emergence of GSs 
in the new bandgaps. The most important result is that in each new bandgap, there is a single family of 
fundamental symmetric GSs which entirely fills the gap. Such GSs near band edges in nonlinear Schrodinger 
equations with periodic potentials were also studied by different approaches in [181 I33| . In Section [5] we 
discuss bifurcations of GSs from a geometric point of view and illustrate this general approach in numerical 
computations for our problem. 

The majority of the newly established GSs are stable solutions of the model equations. In particular, the 
family of fundamental GSs in the new bandgaps is stable even for negative frequencies, when the ordinary 
GSs are unstable [3J [25] . We do not study the stability of GSs in this paper in detail, referring to [55] for 
the verification of the stability by means of direct simulations. 



2. Preliminaries 

2.1. The supergrating eguations 

The standard model of the Bragg grating (BG) in an optical fiber is based on the equations [7] 

iut + iu x + v + (Id 2 + h\u\ 2 ) u = 0, 

(2.1) 

iv t - iv x +u + (\u\ 2 + \\v\ 2 ) v = 0, 
where u and v are amplitudes of the right- and left-propagating waves in the fiber. In these equations, x is 



the coordinate along the fiber, and t is time. A family of exact soliton solutions of (2.1 ) have been found in 

mm- 

Here we will consider a situation where the BG is subjected to periodic modulation in x, with a period 
2n/k, which gives rise to a supergrating. The general form of the extended model reads 

iu t + iu x + [1 — ecos (kx)] v + ucos (kx + 5) u + (Id 2 + ilttl 2 ) u = 0, 

(2.2) 

iv t — iv x + [1 — ecos (kx)] u + /zcos (kx + 5) v + (\u\ 2 + ^\v\ 2 ) v = 0. 

The real parameter e > accounts for the periodic modulation of the BG strength. In real fiber gratings, 
it may be typically 1 cm [51 135j. while the total length of the grating may be up to 1 m. The modulation is 
implemented through the change imposed on the local variation of the refractive index in the fiber's cladding 
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and may also affect the local index felt by each wave. This is accounted for by the perturbation parameter 



fi > and a phase shift 5. In practice, the harmonic superlattice corresponding to Eqs. (2.2) may be 
implemented if the optical interference pattern, which is used to burn the grating into the fiber, is created 
with spatial beatings. We consider stationary GSs, which are found via the ansatz 

u(x,t) = exp(—iut)U(x), v(x,t) = exp(— iuit)V(x). 



Substituting this into Eqs. (2.2 1, we obtain 

ujU + iU' + [1 - ecos(fcz;)] V + fi cos (kx + 5)U+ (\V\ 2 + ±\U\ 2 ) U = 0, 
ujV -iV + [1 -ecos(fcx)] U + ficos{kx + 5)V + (\U\ 2 + ±\V\ 2 ) V = 0, 

where the prime stands for d/dx. 
2.2. Linear analysis 



(2.3) 



We first analyze the bandgap structure of (2.2) by looking for solutions of the linearized system in the 
form of 

u(x, t) = exp (iqx — iujt) U (x), v(x, t) = exp (iqx — iuit) V(x), 



with q being the propagation constant. Thus, we have to find periodic solutions of the linearization of (2.3 1 

(w - q)U + iU' + [1 - e cosikx)] V + u cos (kx + S)U = 0, 

(2.4) 

(lu + q)V -iV'+[l-e cos(kx)} U + ficos (kx + 5) V = 0. 

Note that in the unperturbed problem, with e, /i = 0, the dispersion relation is given by w 2 = q 2 + 1, so that 
we find the well-known gap in the spectrum — 1 < oj < 1. 



Let us consider the case of \i — 0. For e > non-trivial solutions of (2.4) emerge due to parametric 
resonance, similar to the situation in the classical Mathieu equation [37 . Consequently, we expect new gaps 
to open up at points u>± m = ±yl + (mfc) 2 /4, m — ±1, ±2, . . .. We will denote the corresponding gaps by 
1,2,..., while the central gap will be denoted by 0. In the central gap 0, GSs for ^ are found 
as perturbations of the soliton for e, [i = 0, established in [T]. We will use Melnikov's method [T7] to prove 
the existence of gap solitons in Section [3] For the non-central gaps, however, we will apply the averaging 
method to investigate the emergence of (small) GSs for e,fx ^= in Section [4] 

Straightforward perturbation theory for parametrically excited systems yields approximation to the gaps 
for small e > 0. For example, an approximation to the gap with uj close to w±i, i.e. close to the unperturbed 



first gap at e = 0, can be obtained by looking for solutions of (2.4) to lowest order approximation 



U (x) —A + a± cos(kx) + a2 sin(fcx), 
V(x) =A + (3i cos(kx) + 02 sm(kx). 
4 



Using this ansatz we derive at the solvability condition at order 0(e 2 ), 

(c 2 - q 2 - 1) [{u 2 ~q 2 -l~ fc 2 ) 2 - 4fcV] 

= £ 2 [(w 2 - q 2 ) 2 - k 2 (uj 2 + q 2 + 1) - 1] . 

Observe that u> 2 s» (? 2 + 1, since we are interested in the gap close to the unperturbed one with e = 0. 
Moreover, note that at the boundaries of the gaps periodic solutions with wavenumber k or k/2 exist [37 . 
Hence, being interested in solutions with q k/2 we introduce 

P := lu 2 - q 2 - I, Q:=k 2 -Aq 2 , 



and expand (2.5) for small P,Q to obtain 

2fc 2 F 2 + ((fc 2 - 2)e 2 - k 2 Q)P + 2k 2 {q 2 + l)e 2 = 0. 

A new gap appears if this quadratic equation for P has no real solutions. Thus we obtain as the condition 
for the gap, that 

i ■> 

< 2eVk 2 + 4. 



2 -e 2 



k 2 

In terms of the original parameter uj we find that the gaps open up at w±\ = v/l + fc 2 /4 and their widths 
are 

Aoj = = — £. 

w 2 

Similar expressions can be obtained for the higher-order gaps for m > 2 showing that the gap width of 
the band corresponding to m scales with e' m '. We can treat the case of \i ^ similarly. 



These perturbative results can be extended numerically by detecting parameter values, for which Eqs. ( 2.4 ) 
possesses periodic solutions with wavenumber k or k/2. It is well known that these parameters form the 
boundary of the gaps in the spectrum (compare with above) . 

Following [38] we have computed these parameter values as solutions of the corresponding eigenvalue 
problem in Matlab. The results for several combinations of parameter values k and /i are shown in Fig. [T] 
In this figure the shaded and white areas stand for bands and gaps, respectively. It is found that varying k 
merely changes the width of bands and gaps with respect to u>, whereas changing \x leads to more complex 
effects. In particular, for fi =^ gaps opening at e, fi = can vanish again when s or \x are increased. 

3. Existence of gap solitons in the central gap 



In this and next sections we investigate the existence of GSs in Eq. (2.2), using both analytical and 
numerical techniques. These solitons are described by homoclinic solutions to the origin (U, V) — (0, 0) in 
Eqs. Q. 
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Figure 1: The linear spectrum of Eq. l |- 4| in the ((jJ,c) plane. Shaded and white areas are, respectively, bands and gaps. The 
diagrams show the spectrum for a) k = 1, /J, = 0, b) k = 0.5, \x — 0, c) k = 2, /J, = 0, and d) k = 2, fi = 0.5, <5 = 0. 



Before we present the computations let us take a closer look at properties of Eqs. (2.3). It will be 



important that these equations allow for a symmetry reduction by setting V — — U* , where '*' denotes 
complex conjugation. In the following, we will thus only consider the reduced equation 
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iU' +u)U-(\ -ecoskx)U* + ficos{kx + S)U + ^\U\ 2 U = 0. 



(3.1) 



Let U = a + ib, where a,kR. Eq. (3.1 ) is rewritten as 



a' = — [u + 1 — ecos kx + ficos(kx + 5)]b — - (a 2 + b 2 )b, 

3 

b 1 =[ui — 1 + ecoskx + [icos(kx + 5)]a + - (a 2 + b 2 )a. 



(3.2) 



We will consider (3.2) as a dynamical system in the three-dimensional phase space K 2 x S 1 , with S 1 as the 



circle of length 2n/k. 



Note that if j«sin<5 = 0, then Eqs. (3.2) are reversible, that is, invariant under compositions of time 



G 



reversal and the (linear) involutions, 

R\ : (a, b, x) n- (a. —b, —x), R[ 
i?2 : {a, b, x) i — ^ (—a, 6, — x), i? 2 



a,b,x ■ 
a.b.x ■ 



We refer to [2 3) for a general review of reversible systems. Of particular importance for us will be symmetric 
homoclinic orbits, which are mapped to itself under the action of Ri^ or 2 . It is well known that an 
orbit is symmetric if and only if it intersects the invariant plane of the involution. 

Now we assume that < e, /i <C 1 and — 1 < uj < 1, and set u> = cos 9 for some 9 € (0, ir). The cases of 
\u>\ > 1 and fa 1 will be treated in Section [4] and Appendix B respectively. 



For e, (i = Eqs. (3.2 1 become a planar Hamiltonian system, 

a' = -(^ + l)&-^(a 2 + &>, 



b' =(w - l)a + -(a 2 + 6 2 )a, 



(3.3) 



with a Hamiltonian 



i[( w -l) a 2 + (c + l)6 2 ] + 5 (a 2 + fe2)2! 



ff(o,6) 

where a and 6 represent the canonical momentum and coordinates, respectively. The equilibrium at the 



origin is a hyperbolic saddle in (3.3) and has a pair of homoclinic orbits 

cosh(xsin( 



(a±(x),b ± {x)) = ±2 



sin ft cos 



±2\l - sin sin 
3 



2 J cosh(2a: sin 9) + cos 9 ' 
sinh(a;sin#) 



• f 6 \ 

im - 

V2y 



COS ( 



(3.4) 



cosh(2x sin( 

which are symmetric under the involution R\. When e and/or [i are nonzero but sufficiently small, the 



origin is still a hyperbolic saddle in (3.2) and has two-dimensional stable and unstable manifolds which may 
intersect transversely. Such intersection yields transverse homoclinic orbits to the saddle at the origin, which 
persist under variation of the parameters [41] . Along certain curves in the parameter space, however, the 
intersection of the manifolds may become tangential, resulting in bifurcations of the homoclinic orbits. See 
Section [5] for details. 



Here we aim to prove the existence of such transverse homoclinic orbits in Eqs. (3.2) by means of 
the Melnikov method 17 . From formula (4.5.6) in [17], we derive the Melnikov functions M±(xq) for 
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(a±(x),b ± (x)) 



1*00 

M±(xq) =e / [2cj + 3(a±(x) + b±(x))] a±(x) b±(x) cosk(x + xo)dx 



2fi 



a± 



= — e 



-00 
00 



(x) b±(x) cos[fe(x + xq) + 5}dx 



- 2/1 



[2uj + 3(a|. (x) + b±(x))] a±(x) b± (x) sin kx dxj sin kxo 
a± (x) £>± (x) sin kx dx I sin(fexo + S) , 



(3.5) 



where we used the fact that a±(x) and b±(x) are even and odd functions of x, respectively. Substituting 
(3.4) into (3.5) and using the method of residues, we compute 



/oo 
a± (x) b± (x) sin kx dx 
-00 



-QQ 

2 ■ 3 
= - sin J 

3 



sinh(2x sin#) 
[cosh(2x sin 0) + cos 8] 



- sin kx dx 



nk 



cosech 



kn 
2 sin 6 



sinh 



k9 
2 sinO 



and 



[a±(x) 4- b±(x)]a±(x) b±(x) sin kxdx 



- sm° 



sinh(2x sin (9) 
[cosh(2x sin 9) + cos 0] ; 



sin fex dx 



?rfc 



cosech ( 



A'TT 



9 \2sm6 



I k6 \ 
—2 sinh ( — — - ) cos 6 + k cosh 



V 2 sin 



k9 
2 sin 9 



See Appendix A for derivations of (3.6) and (3.7). Hence, the Melnikov functions become 



M±(x ) =- — £ 



2ttA: 



7rfc 2 / /C7T 

cosech - 



\2sin 9 



cosh 



A;0 
2 sin 9 



sin fcxg 



: , / kir \ . / k9 \ . 

-« cosech — ; — - sinh — - — - sinifcxn + d), 
P \2sin0y \2sm0 J v ; ' 



where we used the relation u = cos 9. 
Let us assume that 



,2[i , / k 
¥= 7r or £ # — tanh — 



fe0 



V 2sin 6< 



Then we easily see that M±(xq) has simple zeros at 



7T 27T 
X = Xq, Xq + - mod — , 



where 



xo = arctan 

k 



2/i sin 6 
2/j, cos S + ke coth(fc^/2 sin 9) 



(3.6) 



(3.7) 
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Figure 2: Stable manifold (red) and unstable manifold (blue) of the origin in the Poincare section x = mod 2tt for e = 0.01, 
H = and k = 1: (a) u = -0.5; (b) w = 0.5. 
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Figure 3: Symmetric GSs with V = — U* of in the central gap with k = 1, e = 0.01 and = 0. The solid curve shows the 
real part and the dashed curve the imaginary part of the solution. 



This implies that there exist transverse homoclinic orbits near (a, b) = (a± (x — xq) , b± (x — xq)) and (a± (a; — 
— ir/k),b±(x — xq — ir/k)). In particular, when /xsin<5 = 0, we have xq = so that the first and second 
orbit are symmetric under R\ and i?2, respectively. Thus, we can prove that there exist GSs with V = —U* 



in (2.2| 



To illustrate the above analysis we now introduce the Poincare section {x = mod 2tt} in phase space, 
and consider the return map induced by the flow of ( |3.2[ ). We have computed the stable and unstable 
manifolds of the origin using the continuation tool AUTO 111 with assistance of the HomMap driver [^51 EH] . 
In particular, to compute the stable and unstable manifolds, small segments (~ 10~ 5 ) in the stable and 



unstable subspaces of the origin were taken and solutions of (3.2) with initial conditions on the segments 



were numerically integrated. The same approach was also used for computation of the one-dimensional 
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(b) uj = 0.5 

Figure 4: Asymmetric GSs with V = —U* of f \-.2\ i in the central gap with e = 0.01, k = 1, = 0.01 and 5 = 7r/2. 



unstable manifolds for a three-dimensional Poincare map in [57]. As detected in the theory, the stable and 
unstable manifolds are observed to intersect transversely and GSs could be found. The stable and unstable 
manifolds on the Poincare section x = mod 2tt for [i = 0, u> = ±0.5, e — 0.01 and k = 1 are drawn in 
Fig. [2] Intersections of these manifolds give rise to homoclinic orbits or GSs, and symmetric GSs in (2.2 1 for 
the same parameter values are plotted in Fig. [3| These GSs are not visibly different from the unperturbed 
ones of In Fig. [I] asymmetric GSs in ( |2~2| ) are also plotted for p = 0.01, 6 = ir/2, uj = ±0.5, e = 0.01 

and k = 1. Note that there are different asymmetric GSs which have almost the same shapes as those in 
Fig.|4]but have peaks at a position shifting by n/k. Plots of GSs throughout this paper show the real part 
of the solutions as a solid curve, whereas the imaginary part is plotted as a dashed curve. 



4. Existence of gap solitons in non-central gaps 



In order to investigate the existence of GSs of (2.2) with V = —U* in the non-central gaps, we shall use 
the higher-order averaging method [2S][2!|]. Let us describe a general setup for the averaging analysis. 



We seek homoclinic solutions to the origin of (3.2 ) near uj = uj± m in the limit of small e, fi. Note that for 



e = fi = the origin is an equilibrium in (3.2), at which the Jacobian matrix has a pair of purely imaginary 
eigenvalues A = ±i\/uj 2 — 1 when uj s» oj± m . Let j and I be non- negative integers which will be determined 
for each gap later. Define a detuning parameter Q by 

e 3+1 il = uj — uj± m . 



Let 



/ 



$(ar) 



mk 

cos 



in 



k 



. mk 

, sin x 

\2(cu ±m + l) 2 
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2(w±r. 



mk 

mk 

cos x 

2 



mk \ 
~2 



which is the fundamental matrix of the linearization at the origin of the unperturbed system (3.3) with 
lj = io± m . Using the transformation 

= e (»+l)/2 $(a .) 



k 
m 



\ 



\2(w ±m + l) ? V 



(4.1) 



in (3.2 1, one obtains 




1), 



■e>VL + coskx - [icos(kx + S) - |e'(a 2 + b 2 )]b ' 



^ [e j n + cos kx + (i cos(kx + S) + \e l (a 2 + b 2 )]a 



(4.2) 



where /i — fi/e and a and 6 are represented by £ and rj via (4.1 1. 



Homoclinic solutions in (4.2), equivalently in (3.2), can be approximately obtained by the averaging 
method. These homoclinic solutions correspond to GPs of (2.2) with V — —U* in the non-central gaps. 
While in the first gaps the standard (i.e., first-order) averaging method is sufficient, we have to use 
the higher-order averaging method for the higher-order gaps. This difficulty occurs because the potential 
considered here consists of the pure harmonic component coskx. But it is important to note that higher- 
order averaging is also needed for more general potentials with higher Fourier components, if these higher 
components are small compared to the first one. This is the reason for us to present the higher-order method 
here. 

The necessary lengthy computations can be easily performed with computer algebra systems. Specifically, 
we have used the program haverage .m for the computer software Mathematica |42j . developed in [461148] (see 
also [35]). We exemplarily present the results for the first three gaps 1 , 2^ and 3 ± . Similar second-order 



averaging analyses can also calculate the boundaries of the central gap (see Appendix B ) 



4-1. Averaging in the gaps 1^ 



For the computations in the first gap at uj±i = ±yl + fc 2 /4, we set m — 1 and j = I = 0. The standard 



averaging procedure for (4.2) yields the first-order averaged system 

fLsinS w±i(l — 20) + ficosS 



n 



k ^ ' k 

Co ±i (1 + 2il) + fl cos S /2sini5 



k 



?7-7±i(£ 2 +V 2 )V 
?/ + 7±i(£ 2 + ?7 2 )£ 



(4.3) 



where 



7±1 = 2fc3 (w±1 + l)3> - 
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Eqs. (4.3) are easily obtained by averaging the right hand side of (4.2 1 with respect to x over [0, 2tt/A;] and 



such computations can be done by the Mathematica program haverage.m. In the averaged system (4.3 1 
with /2cos<5 + uj±i 7^ 0, we perform a rotational transformation 




to obtain a system of a more convenient form 

2w±i 



(fti±-fi)r?-7±i(f + 



n' = e 



k 
k 



(f) 1± + fiK + 7±i(C 2 + '7 2 K 



where 



and 



n 1± = 



U!±l 



2 | — ) cos S + 1 



for oj±i(/j,cos5 + uj±i) > 0; 



7i + — for w±i(/_tcos(5 + w±i) < 



with 



1 / ^sin<5 

7i = — arctan — 

2 \ cos o + uj± i 



€ [— 7T,7r]. 



(4.4) 



Note that Eq. (4.4) is still valid for (icos5 + uj±i = if 9 = ir/A is taken when uj±ifxsmS > and if ( 



-tt/4 



is taken when w± i /i sin <5 < 0. The system (4.4) is Hamiltonian with a Hamilton function 

' CJ±i 



(Oi± + n)i 2 - (n 1± - ] + i 7± i(l 2 + ?? 2 ) 2 



for which the origin (£, fj) — (0,0) is a saddle if fl e (— fii±,f2i±), i.e. 



w e (w±i - efli±, lj±i + eOi±) 



(4.5) 



We easily see that if 



(i = e\ 1 



8 = (resp. S = n) 



(4.6) 



then = (resp. f2i_ = 0) and a closing of the gap region 1 + (resp. 1 ) occurs. See also Fig. [TJd) . If 
5 0, 7r, then J7i± ^ 0, and thus such a closing does not occur. 



Let us assume that fii-t ^ 0, i.e., condition (4.6) does not hold. A straightforward analysis of the level 



sets of Hi shows the existence of two homoclinic orbits in (4.4) and hence in (4.3). Moreover, we can obtain 
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analytical expressions of the homoclinic orbits as 




2w±ifti± . . 1 
sin © sin 



sinh[(2w±if2i± sin <p/k)ex\ 



3kj. 



±1 



2w±if2i± . 

— — sm <p cos 

3fc7±i 



2 cosh[(4cj ; ti0i± sin <p/k)ex] + cos (/> ' 

1 cosh[(2o;±ini± sin^/fc)ea;] 

2 cosh[(4a;±irii± sin<^)/A;)ea:] + cos0 



(4.7) 



where f2/fix± = cos (p. 

It is easy to see that, since the coefficient in front of ij 2 in Hi is negative, the homoclinic orbits intersect 
the ry-axis and are symmetric about it, forming a figure-of-eight configuration. In general, it is very difficult 



to make precise statements about the symmetry of homoclinic orbits in (4.3 1. For some special cases, 
however, we have the following results: 



(i) When ft > and 6 = 0, the homoclinic orbits in the averaged system (4.3) are symmetric about the 
77-axis if p, < — u}±\ and about the £-axis if p > —u>±i in the negative gap, while they are always 
symmetric about the ry-axis in the positive gap. 



(ii) When [i > and 5 = tt, the homoclinic orbits in (4.3 1 are symmetric about the ?7-axis if p < ui±i and 
about the £-axis if p > u±i in the positive gap, while they are always symmetric about the 77-axis in 
the negative gap. 



(iii) When /i = 0, the homoclinic orbits in (4.3) are symmetric about the 7/-axis in both the positive and 
negative gaps. 



From the existence of homoclinic orbits in the averaged system (4.3 ), we claim that there exist homoclinic 



orbits in (3.1) and hence stationary GSs in (2.2). For this consider (4.2) and (4.3) as dynamical systems 



defined in the phase space R 2 x S 1 . Recall that we set m = 1 and j = I = in (4.2). We introduce a cross 



section S at x — mod 2n/k and define the Poincare maps ip, ip : S —> S for (4.2) and (4.3), respectively. 
For both t/j and %j) the origin (£, 77) = (0, 0) is a fixed point and has stable and unstable manifolds. Moreover, 



the stable and unstable manifolds, W s (0) and W u (0), coincide along the homoclinic orbits in (4.3) for ip. 
By the averaging theorem (see Theorem 4.1.1 of [T7]), the stable and unstable manifolds, 1F S (0) and W u (0), 



for ip are 0(e)-close to W s (0) and ^"(0). Since Eq. (4.2) as well as Eq. (4.3) are Hamiltonian, the Poincare 
maps ip, ip are area-preserving. 

Now assume that the stable and unstable manifolds VK S (0) and W u (0) do not intersect, and go inward 
and outward, respectively, as shown in Fig. [5j We want to show that this yields a contradiction and define 
a region S, which is encircled by two parts of W s (0) and W u (0) and a line connecting two points pi and P2 
on VF s (0) and W u (0). The region S is mapped by the Poincare map ip to another region S' = ip(S) which is 
encircled by two different parts of VF s (0) and W""(0) and a curve connecting two points ip(pi) and ip(j>2) on 
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Figure 5: Stable and unstable manifolds W s ' u (0) (solid lines) which are close to W s ' u (0) (broken lines) but do not intersect. 
The regions S and S'\S with S' = ip(S) are hatched distinctly. 



W /s (0) and ^"(0). It is obvious that S' is larger than S and contradicts the fact that ip is area-preserving. 
Similarly, we cannot assume that W s (0) and W U (Q) go outward and inward, respectively. Therefore, we 



can prove that W s (0) and W u (0) intersect and there exists a homoclinic orbit in the full system (4.2) and 



hence in (3.1). On the other hand, Eq. (4.2) can be averaged up to 0(e n ) for any integer n > 1 but has 



very rapid oscillations compared with the unperturbed homoclinic orbits (4.7) which vary very slowly with 
0(l/e). Applying a result by Neishtadt [3D] (see also [ID]), we see that the splitting distance between W s (0) 
and W"(0) is exponentially small with respect to e at most. Thus, there exist stationary GS solutions to 



the system of equations (2.2) in the gap regions (4.5). 



Moreover, assume that /isin<5 = 0, i.e., fj, — or S = 0, n, so that Eqs. (3.2) is reversible with respect to 
Rj and R'j, j = 1,2. Then it follows from the above arguments that the homoclinic orbit in the averaged 



system (4.3) is symmetric about the £-axis (resp. ry-axis). By the persistence of symmetric orbits in reversible 



systems, there exist symmetric homoclinic orbits about the £-axis (resp. ?7-axis) in the full system (4.2) and 



with respect to R\ and R[ (resp. i?2 and R' 2 ) in (3.2) and hence in (3.1). In particular, the homoclinic 



orbits in (3.1) have an even (resp. odd) real part a and odd (resp. even) imaginary part b. 



We can use similar arguments to the central gap when condition (3.8) does not hold, and prove that 
there must be homoclinic orbits. However, the splitting distance between the stable and unstable manifolds 
is not exponentially small with respect to s but 0(e 2 ) at most since both of the perturbations in (3.2) and 



unperturbed homoclinic orbits (3.4) vary with O(l) 



In addition, we can detect a non-transverse homoclinic orbit in (3.2) when 5 — n and 

2/i 



k 



tanh 



k9 
2 sinfl 



(4.8) 



as follows. Fix S = it. Then the Melnikov functions M±(xq) have graphs as shown in Fig. [6j depending on 
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M t (x ) 



8<8 




8=8 o 




8>8 













Figure 6: Graphs of the Melnikov functions M±(xq) for S = tt. 




Figure 7: Stable and unstable manifolds W s ' u (0) for the associated Poincare map: (a) e < so; (b) £ > £o\ ( c ) £ = £0- 



£ < e 0l £ = £ o or £ > £ o- Since M±(x ) represent signed measures of their distance (see, e.g., section 4.5 



of [H]), noting that b±(xo) = and Eqs. (3.2) is reversible, we can draw the behavior of the stable and 
unstable manifolds W s ' u (0) for the Poincare map of (3.2 1 as shown in Figs.|7ja) and (b). Using the fact that 



the Poincare map is area-preserving and applying a discussion similar to the above one for (4.2), we see that 
W s (0) and W u (0) coincide for some e near e , say Eq, as shown in Fig.[7]jc). Thus, the point (6,e) = (tt, £q) 
is very degenerate in the parameter space. This also suggests that new interesting behavior may occur at 



the point in the original PDEs (2.2) 



In Fig. [8] we show the results of numerical computations of the stable and unstable manifolds by AUT097 
with HomMap near the homoclinic bifurcation point e = Eq ps 0.010806 for /i = 0.01, u> — 0.5 and k = 1. 
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-0.8 -0.4 0.4 0.8 

a 



Figure 8: Stable and unstable manifolds of the origin on the Poincare section x = mod 2-k when fi = 0.01, ui = 0.5, S = n 
and k = 1: (a) e = 0.0088063; (b) e = 0.0108063; (c) e = 0.0128063. The red and blue curves represent the stable and unstable 
manifolds, respectively. 




-1 -0.5 0.5 1 

CO 

Figure 9: Homoclinic bifurcation curves for Eqs. with k = 0.5, 1,2. The parameter values (oj,e) of Fig. [8jb) with k = 1 

is plotted as •. 
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We see that as predicted by the theory, the invariant manifolds almost coincide at e ~ Sq while they split 
when e is different from bq. Thus, the theoretical prediction for the homoclinic bifurcations is very precise. 
Figure [9] also shows the approximate homoclinic bifurcation curves given by Eq. (4.8), on which the stable 
and unstable manifolds coincide, in the (to, e/^)-space for k = 0.5, 1, 2. 

4-2. Averaging in the gaps 2 ± 

We set ui — 2 and j = I = 1, and perform a second-order averaging method to prove the existence of gap 



solitons of (2.3) with V = —U* in the second gaps appearing at uj±2 = ±vl + k 2 . To avoid the complexity 



of expressions, we give our results only for fi = in the gaps of higher orders. The second-order averaged 
system has been obtained using the Mathematica program haverage . m, and it reads 

^2 



t . 2 {a±2 [(6(w±a + m - 5)fc 2 + (w ±2 + l) 2 
+ P±2 [k 2 e + (w ±a + 1) V] } V, 



(4.9) 



— ^ {a ±2 [(6(w± 2 + m + 1)^ 2 - 5(^ ±2 + l) 2 ] 

+ P±2 [k 2 e + (U±2 + 1) V] } e, 



with 



a±2 



±2 



fc 2 + (LU ± 2 + l) 2 

" 12fc 2 (w ±2 + 1) ' 

3[3fc 4 + 2(u ±2 + l) 2 fc 2 + 3{uj ±2 + l) 4 ] 
16fc 2 (a; ± 2 + l) 2 



Note that a+2 = — a-2 > and /3±2 > 0. Eq. (|4.9|) is a Hamiltonian system with the Hamilton function 
H 2 (LV) =e 21 



a± 1/1 {[6(w ±1/1 + 1)0 + l]fc 2 - 5(^1/1 + l) 2 } 
2( W±1/1 + 1)2 



+ 



«±x/i{[6(w±i/i + 1)0 - 5]fc 2 + (w ±1/1 + l) 2 } 2 



2fc 2 



+ ,, 2 / ±l \ 1 ,2 [ fc2 e 2 + (^i/i + i: 



V] 2 



Let 



AO2 



AO 2 . 



2k 2 - y/k 2 + 1 - 1 2fc 2 - 5Vfc 2 + 1 + 5 



3A; 2 (yFTT+i) 3/c 2 (VFTT-l) ' 
2fe 2 + 5Vfc 2 + 1 + 5 2fc 2 + Vk 2 + 1 - 1 



3fc 2 (7F 



1 + 1) 



3fc 2 (7fc 2 



1-1 



The origin is a saddle of (4.9 1 and has a pair of homoclinic orbits if O 6 (A0 2 , AO 2 ) and O e (— AOf,, — A0 2 ) 



when the signs '+' and ' — ' are taken in the subscript of uj in (4.9), respectively, i.e 



u G + e 2 A02, uj+2 + ^AO 2 .) 
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or 

Again, we find the averaged system to be Hamiltonian with a saddle equilibrium at the origin, and a 
study of the corresponding level set of the Hamiltonian shows the existence of a figure-of-eight of homoclinic 
orbits. These orbits intersect the 77-axis. By the same argument as above, this implies the existence of gap 



solitons for (3.1 ) in the gaps 1 having an even a-component and an odd ^-component. 



4-3. Averaging in the gaps 3 4 



We set to = 3, j = 1 and I = 2, and perform the third-order averaging procedure for (4.2) using the 
Mathematica program haverage . m to obtain the averaged system 



-2 



where 



= - { t 32 ( w ±3 + ^ - 9]fc 2 - 4( W±3 + l) 2 } n 

+ ^5 {«±3 ~ /3 ±3 [9fc 2 £ 2 + 4( w±3 + 1) V] } V, 
V' = 4(cj J" +1)2 P±3 {[32( W±3 + 1)0 - 9]fc 2 - 4(c ±3 + l) 2 } £ 
+ 4( ^J 3 + 1)2 {"±3 + /3±s [M 2 e + 4( w±3 + 1) V] } 



= 9fc 2 +4(^ ±3 + l) 2 
P±3 ~ 64fc 2 ( w±3 + l) ' 

_3[9fc 2 + 4(w ±3 + l) 2 ][81/c 4 - 56(cj ± 3 + l) 2 fc 2 + 16(cj ±3 + l) 4 
tt±3 ~ 512fc 4 (w ±3 + l) 2 

81fc 4 + 24(cj ±3 + l) 2 fc 2 + 16(w± 3 + l) 4 



^ 3 64fc 2 (w ±3 + l) 2 

Note that p + z = — ps > and a±3,/3± 3 > since 

81fc 4 - 56(w ±3 + l) 2 k 2 + 16(w ±3 + l) 4 
= [9fc 2 - 4(w ±3 + l) 2 ] 2 + 18(w ±3 + l) 2 > 0. 

Let 

Q = ± 3 + ev, 

where 

= 9fc 2 +4(^ ±3 + l) 2 
±3 32( W±3 + l)fc 2 ' 



Note that Q +3 > and < 0. Eq. (4.10) becomes 



e = - ^{(7±3^ - a ±3 ) + P± 3 {9k 2 C 2 + 4( W±3 + 1)V]K 

rf = 4(cj J 3 +1)2 {(7±3^ + a± 3 ) + /3 ±3 [9fc 2 ? 2 + 4( W±3 + 1) V]R, 
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(4.10) 



(4.11) 



where 



7±3 = 32p± 3 (u±s + l)fc 2 = -[9k 2 +4(u ±3 + l) 2 ] > 0. 



Eq. (4.11 1 is a Hamiltonian system with the Hamilton function 



H^,'n)=E a \— — -£ + — — r\ 



(W±3 + I) 2 
/?±3 



144fc 2 (w ±2/1 + 1) 



18fc 2 

[9fc 2 e 2 +4( w±2/1 + l)V] 2 



(4.12) 



We easily see that 



o+3 = a_3 = 3(fc 4 + 4) 

7+3 7-3 



16fc 4 



and set 



A^3 = 



3(fc 4 + 4) 
16fc 4 



The origin is a saddle and has a pair of homoclinic orbits if v £ (— AS1 3 , Ails), i- e -> 

cj e (cj ±3 ± e 2 fi 3 - e 3 An 3 ,uj ±3 ± e 2 fl 3 + e 3 Aft 3 ) . 



Especially, in the regions (4.3) the coefficients of £ and 77 in the Hamiltonian (4.121 are positive and 



negative, respectively, so that the pair of homoclinic orbits draws a vertical figure-eight' in the (£,?7)-phasc 
plane. 

Similar computations reveal the existence of gap solitons in the higher order gaps. These solutions are 
odd in a (and even in b) or vice versa, depending on whether the gap number is odd or even. A general 
proof of this is beyond the scope of this paper. In the following we illustrate and extend this perturbative 
results by numerical computations. 

4-4- Numerical studies 



We computed the stable and unstable manifolds of the origin and homoclinic orbits in (3.2), i.e., GSs with 



V = —U* in (2.2), using the computer package, AUTO with HomMap. In our computations we found GSs due 



to the transverse intersection of the stable and unstable manifolds, thus confirming the analytical results. 
For several parameter values, the stable and unstable manifolds for the Poincare section x — mod 2n 



are shown in Fig. 10 and GSs in (2.2) are plotted in Figs. 111113 In addition to the solutions of (3.2), the 



figures include homoclinic solutions of the averaged systems (4.3), (4.9) and (4.11), rescaled according to 



the transformation (4.1). As can be seen, they provide a good match to the envelope of the solutions to 



(3.2). In all figures presented, the parameter k has been set to k = 1. 

In closing this section, we return to the parameter values of Fig. [ljd) . In this case a new bandgap opens 



up for small e, but closes at e = l/y/8. See also formula (4.6). Our computations reveal that GSs exist in 



the gaps both above and below this parameter value. In Fig. 14 we show solution profiles at e = 0.2 and 
e = 0.8. 
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Figure 10: Stable and unstable manifolds of the origin on the Poincare section x = mod 2ir when fi = and k = 1: (a) 
a; = -1.118, £ = 0.045; (b) ui = 1.118, e = 0.045; (c) uj = -1.42, e = 0.4; (d) u = 1.42, e = 0.4; (e) u = -1.9, s = 0.6; (f) 
w = 1.9, e = 0.6. The red and blue curves represent the stable and unstable manifolds, respectively. 

5. Bifurcation of solitons 

We further discuss the behavior of GSs under variation of the system parameters. When the parameter 
values are varied, homoclinic tangencies, i.e., homoclinic bifurcations, can occur, leading to the creation of 
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(a) w = -1.118, £ = 0.045, fj, = 
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(b) w = 1.118, e = 0.045, n = 
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-0.40- 




-15 15 

(c) u = -1.118, £ = 0.4, = 



-15 15 

(d) uj = 1.118, £ = 0.4, fj, = 



Figure 11: GSs with V = — U* of j2.3| in the gaps 1^, 2^ and for k = 1, symmetric with respect to i?i or i?2- Gray curves 
in the panels correspond to solutions of the associated averaged system, |4.3| , |4.9[ l or | |4.11[ l. 

additional GSs and bifurcations of GSs. We have found such bifurcations in all gap regions including the 
central one. 



As an example, we plot bifurcation diagrams of GSs for Eqs. (2.3 ) in the gaps 1 and when ui = —1.118 



or 0.5 and 5 = n/2, in Fig. 15 Here (a(0),6(0)) represents the point at which each GS crosses the section 
S = {x = mod 2ir}. In the diagrams we observe saddle-node and pitchfork bifurcations of GSs. The 



occurrence of these bifurcations is explained from the behavior of the Poincare map of (3.2) for the section 
E, as follows. 

Let us first assume that /i ^ 0, such that the Poincare map is not reversible. In this case, a tangency 



between the stable and unstable manifolds as shown in Fig. 16 '&) leads to a saddle-node bifurcation, in 
which two new GSs are created. See Fig. |15[ c) for a numerical evidence. 

Now suppose that fi = and let us choose a cross-section, such that the Poincare map is reversible with 
respect to R\ and i?2- We can distinguish two different bifurcation scenarios. 
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(g) to = -1.9, e = 0.6, fi = (h) w = 1.9, e = 0.6, \i = 

Figure 11: Continued. 




-80+7C 



(a) cj = -1.118, e = 0.045, /x = 




-80+jc 



(b) w = 1.118, e = 0.045, /i = 



Figure 12: GSs with V = — U* of |2.3i in the gaps 1^ for k = 1, symmetric with respect to or H^. Gray curves in the 
panels correspond to solutions of the averaged system | |4.3| . 



If the stable and unstable manifolds have a cubic tangency on the invariant plane of R\ or i?2 , but intersect 
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(a) uj = -1.118, e = 0.045, fi = 0.045, 6 = tt/2 (b) u = 1.118, e = 0.045, /u = 0.045, S = tt/2 



Figure 13: Asymmetric GSs with V = —U* of j2,3[ ) in the gaps 1^ for k = 1. Gray curves in the panels correspond to solutions 
of the averaged system i|4.:!|. 




(a) s = 0.2 



(b) e = 0.8 



Figure 14: Symmetric GSs with V = —U* of \2.3\ in the gaps for w = —1.4, fi = 0.5, 8 = and k = 2. 



this plane transversely as shown in Fig. 16 'b), then a pitchfork bifurcation of homoclinic orbits takes place. 



When the parameters are varied, a pair of homoclinic orbits appear, while the original homoclinic orbit still 



exists due to the persistence of symmetric orbits. See Figs. 15 a) and (b). The original GS is symmetric, 
but the two new orbits are asymmetric. 

On the other hand, if the stable and unstable manifolds are tangent to the invariant plane at their 



quadratic tangency as shown in Fig. 16 ^c), then a saddle-node bifurcation of symmetric GSs occurs. See 



Fig. 15 ^b). Note that a saddle-node bifurcation of asymmetric GSs can occur even though fi = 0. 



We illustrate these general ideas with computations for equation (3.2). Figure 17 shows the stable and 



unstable manifolds of the Poincare map of (3.2 ) on the section S for the parameter values when saddle- node 



or pitchfork bifurcations occur in Figs. 15 
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Figure 17 'a) corresponds to Fig. 16 'a) and depicts the situation when the map is not reversible (/i > 0) 
and the stable and unstable manifolds have a quadratic tangency. Similarly, in Figs. |17[ b) and (c) the 
Poincare map is reversible, and the manifolds have cubic and quadratic tangencies on the invariant plane 
b = as in the sketches in Figs. 16 'b) and (c). 



Figures 18 and 19 show the profiles of GSs born at the bifurcations in Figs.|l5[a)-(c) and (d), respectively. 



In particular, the GSs in Fig. 19 are bound states with two humps. Standard results about transverse 
homoclinic orbits explain the existence of multi-pulse GSs, see |41j . and we also find such GSs in the other 
gap regions. 

We note that bifurcations of symmetric homoclinic orbits in reversible systems that are autonomous have 
been discussed in [6) 122] . Our analysis of the Poincare map also yields comparable bifurcation results in 



systems with periodic forcing like Eqs. (3.2). 
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Figure 16: Homoclinic bifurcations in the Poincare map. Panel (a) depicts the situation that leads to a saddle-node bifurcation 
of asymmetric GSs. Panels (b) and (c), respectively, explain mechanisms for pitchfork and saddle- node bifurcations of symmetric 
GSs. 

6. Conclusions 

In this paper we have studied the existence of GSs (gap solitons) in a model of the BG (Bragg grating) 
subjected to periodic modulation along the fiber, thus creating the supergrating. The corresponding math- 
ematical model amounts to a parametrically forced extension of the standard BG model with the additional 
symmetry-breaking forcing term. 

As was known previously [5], the supergrating opens up a set of new bandgaps. Using the Melnikov 
method and averaging techniques, we have demonstrated that families of GSs exist robustly in the new 
gaps. We have also analyzed bifurcations of GSs and showed that, in addition to the fundamental GSs, the 
model also supports families of bound states of the solitons. 

Many of the GSs found in the model are completely stable, as was verified by means of direct simulations 
for some of them [49] . The present paper extends the results of [49] and gives detailed mathematical proofs 
for them, although we did not analyze the stability problem here in detail. In particular, we have discussed 
the effects of the symmetry-breaking forcing term and demonstrated that it can lead to the closing and 
re-opening of the bandgaps. In the presence of the new term, GSs exist robustly (as generic solutions) in all 
the bandgaps, but, as the symmetry of the equation is broken, they undergo many bifurcations. 
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Figure 17: Stable and unstable manifolds of the origin on the Poincare section x = mod 2tt when k = 1 and 8 = n/2: (a) 



-1.118, e = 0.3, m = 0.0326551; (b) w = -1.118, e = 0.234076, ^ = 0; (c) ui = 0.5, e = 3.35338 x 10" 



0. The red 



and blue curves represent the stable and unstable manifolds, respectively. In plate (a) the invariant manifolds have quadratic 
tangencies inside the circles. 
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Figure 18: GSs with V = -U* of ( |2.3| born at the bifurcations of Fig. |l5|»-(c) for e = 0.4, u! = -1.118, k = 1 and fj, = 0. 
The GSs in plates (a) and (b) are symmetric and born at the saddle-node bifurcation while the GSs in plate (c) is asymmetric 
and born at the pitchfork bifurcation. 
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Figure 19: Symmetric, two-pulse GSs with V = 
uj = 0.5, pu = and k = 1. 
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Figure A.l: The rectangle contour F. 
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Appendix A. Derivations of equations (3.6) and (3.7) 

Let 



sinh(2a; sin 9) 



[cosh(2x sin 6) + cos 0p 



and let 



fj(x) sin kxdx 



(A.l) 



for j = 2, 3. To estimate (|A.l|), we consider the complex integral 

fj(z)e ikz dz, 



27 



where T is a rectangle contour, —R < x < R and < y < ir/s'me, depicted in Fig. A.l with R > 0. The 
complex function fj(z) is singular only at 

7T±(9 



2 sin 6> ' 

which are jth-order poles, inside of T. By the residue theorem, we have 



f j (z)e ik *dz = 2iri (p^+P^), 
where p^. is the residue of fj(z) sinfcz at z = i(^±. We compute 



(±) ke k <± (±) 

P2 =T A . 3- , P 3 



4 sin J (9 ' 



so that 



and 



f 2 {z) e dz = i ■ 3 exp 



16 sin 

kir \ ainh f kd 
sin 3 6 ""^ V 2sin6» J Sm ySsinfl 



7rA- 



7tA' 



h{z) e lKZ dz =z— ^-exp 
4 sin 



—2 sinh 



kix 
2 sin i 



2 sin 6 



cos 6> + k cosh 



fc0 
2sin( 



On the other hand, we can write 



fj(z)e ik * dz 



fj( X y kx dx + i 



7T j sin 9 



f j (R + iy)e ik ( R + i ^dy 



+ / fJx + in/ sm6)e lk{x+l7r/ sin9 Ux 



f f j (-R + iy)e 

J 7T / sin 



ik{ — R-\-iy) 



dy. 



Since by fj(x + in/ sin 6) = fj{x) 

-R f kir \ r R 

f 3 (x + in/sm6y Hx+l7r/sinS) dx = -exp[-—) / fj(x)e ikx 

V smOJ J_ R 



dx 



and the second and fourth integrals in (A.3) tend to zero as R — ¥ oo, we have 

2m (p { - ] + p[ +] 



f 3 {xy kx dx= f ^ 



where we used (A.2|. Hence, we obtain 

h = 

and 



1 — exp I — ; 



cosech ( — ; — - ) sinh ( — '- 



I3 = — W \ cosech 
8 sin 5 (9 



2 sin J 6 

kit 
2 sinO 



ke 



V 2 sin 9 J V 2 sin t 



-2 sinh ( 



ke 

V 2 sin0 



cos e + k cosh 



ke 

2sin( 



from which Eqs. (3.6) and (3.7) follow. 



(A.2) 



(A.3) 
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Appendix B. Calculation of the boundaries of the central gap via averaging 

Applying the third-order averaging method, we calculate the boundaries of the central gap. Again, the 
Mathematica program "haverage.m" is used to obtain the third-order averaged systems. 



Appendix B.l. Right boundary 



Let oj = 1 — e 2 v and pL — efi. Using the transformation (a,b) — e 3 / 2 7/), we rewrite (3.2 1 as 

f = - 2e 1/2 r/ + e 3/2 [coskx - flcos{kx + 8)]r) + 0(e 5/2 ), 

rf =e 1/2 [coskx + fLcos(kx + 5)}£ + e 3/2 (- v £+fifA +0{e 5/2 ). 



(B.l) 



The third-order averaged system for (B.l I is obtained as 



? = -2eV\ rf 



-3/2 



7^(1 + M 2 + 2£ COS <J)W^ 3 



(B.2) 



where e 1 / 2 was taken as the small parameter. In the averaged system (B.2|, the origin is a saddle if 



v > —(1 + pi 2 + 2/icos<5). 

Hence, we have the approximate right boundary of the central gap, 

6e 



u) = 1 —5- (1 + u + 2/^ cos 5). 

k 



Appendix B.2. Left boundary 



Let lu = — 1 + e 2 v and /Lt = e/i. Using the transformation (a, b) = (e 3 ' 2 £ , erf), we rewrite (3.2 1 as 



£' =e 1,2 [coskx -pcos(kx + 5)}n - e 3/2 \ vr) + -rfj + 0(e 5/2 ), 
rf = - 2e 1/2 £ + e 3/2 [cos kx + p, cos(kx + 5)]£ + 0(e 5/2 ), 



(B.3) 



The third-order averaged system for (B.3) is obtained as 



e = e 3 ' 2 



- ( " - p(l +M 2 - 2/xcos(5) U - ^if 



i = -2£ 1/2 £, 



in which the origin is a saddle if 



v > 7-7(1 + (J 2 — 2/jcos<5). 



(B.4) 



Thus, we have the approximate left boundary of the central gap, 

6e 2 



lo = 1 



(1 + y 2 - 2/icosS). 
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